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ABSTRACT 

In this paper, we study the distribution functions that arise naturally during self- 
similar radial infall of collisionless matter. Such matter may be thought of either 
as stars or as dark matter particles. If a rigorous steady state is assumed, then the 
system is infinite and is described by a universal distribution function given the self- 
similar index. The steady logarithmic potential case is exceptional and yields the 
familiar Gaussian for an infinite system with an inverse-square density profile. We 
show subsequently that for time-dependent radial self-similar infall, the logarithmic 
case is accurately described by the Fridmann and Polyachenko distribution function. 
The system in this case is finite but growing. We are able to embed a central mass 
in the universal steady distribution only by iteration, except in the case of massless 
particles. The iteration yields logarithmic corrections to the massless particle case 
and requires a 'renormalization' of the central mass. A central spherical mass may be 
accurately embedded in the Fridmann and Polyachenko growing distribution however. 
Some speculation is given concerning the importance of radial collisionless infall in 
actual galaxy formation. 

Key words: Cosmology:theory - Dark Matter - large-scale structure of Universe - 
Galaxiesdormation - galaxies:haloes - galaxies: bulges - gravitation 



1 INTRODUCTION 

The relation between the formation of black holes and 
of galaxies has developed into a key astrophysical ques- 
tion, from the early papers by Kormendy and Richstone 
( |Kormendy fc Richstone 1995[ ), to the more recent discov- 
eries by Magorrian et al. (Magorrian et al. 1998[), Fer- 



rarese and Merritt (Ferrasc fc Merritt 2000[) . and Gebhardt 
et al.QGebhardt et al., 2000| ). Recent papers ( Graham 2004 ; 
|Kormendy fc Bender 2009[ ) also establish a correlation be- 
tween black hole mass and bulge luminosity deficit. These 
papers establish a strong correlation between what is es- 
sentially the black hole mass and the surrounding stellar 
bulge mass (or velocity dispersion). The origin of this pro- 
portionality, which extends well beyond the gravitational 
dominance of the black hole, remains uncertain. But it is 
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generally taken to imply a coeval growth of the black hole 
and bulge. 

Various proposals have been offered to explain the black 
hole mass-bulge mass proportionality as a consequence of 
the AGN phase. There is as yet no generally accepted sce- 
nario although a kind of 'auto-levitation', based on radia- 
tive feed-back from the accreting black hole to the star- 
forming gas that in turns limit accretion, is plausible. In any 
event there remains the question of the origin of the black 
hole seed masses. In some galaxies at very high red shift 
the inferred black hole masses are already of order 10 9 M@ 
(Kur k~et al., 2007[ e.g.) after about one Ga of cosmic time. 
This may require frequent, extremely luminous early events 
d Walter et al., 2009[ e.g.), or it may suggest an alternate 
growth mechanism. 

The latter possibility is reinforced by the detection of 
a change in the normalization of the black hole mass-bulge 
mass proportionality in the sense of relatively larger black 
holes at high red shift (Maiolino ct al. 2007, e.g.). As sug- 
gested in that paper it seems that the black holes may grow 
first, independently of the bulge. The collisionless matter 
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that we invoke might be stars or it might be the dark mat- 
ter itself. 

Recently (Peirani & do Freitas Pacheco 2008) have 
studied the possible size of the dark matter component 
in black hole masses. By assuming that the dimensional 
or 'pseudo- phase-space density' (H2006] i.e. Henriksen 
2006b, e.g.) is strictly constant they convert the relativistic 
accretion of the dark matter into an adiabatic Bondi flow 
problem and obtain the resulting accretion rate. Then by 
adopting the mass proportionality between bulge and black 
hole and fitting boundary conditions from cosmological 
halo simulations, they deduce that between 1% and 10% of 
the black hole mass could be due to dark matter. 

If we accept this result at face value, a seed mass of 
say 1O 6 M0 could have grown from dark matter. It would 
now be part of a super massive black hole that subsequently 
grew in the AGN phase. Some seeds may be primordial. As 
early as 1978, through fully general relativistic numerical 
collapse calculations, (Bickncll & Henriksen 1979) predicted 
primordial black hole masses in the range 10 2 to 1O 6 M . 

Stellar density cusps surrounding black holes have 
been studied extensively previously. Classic studies by 
i|Peebles 1972[) and by (|Bahcall fc Wolf 1976|) dealt with 
the problem of feeding the black hole from a filled loss 
cone (nearly radial orbits). In addition to these diffusion 
studies, Young (Young 1980) explored the cusps produced 
by the adiabatic growth of a black hole in a pre-existing 
isothermal stellar environment. This was extended by 
( IQuinlan et al 1995[ ) and by (jMacMillan fc Henriksen 2002 P 
to more general environments. In a cosmological con- 
text, Bertschinger ( [1985] ) also studied the growth of a 
central black hole by radial infall. The conclusions were 
that the black hole induced cusps were never flatter than 
r -1 ' 5 (the isothermal and cosmological case) and that no 
black hole mass-bulge mass correlation was established 
(MacMillan & Henriksen 2002). The latter conclusion has 
spurred the investigation of coeval dynamical growth of 
the black hole and bulge in contrast to adiabatic growth 
(MacMillan & Henriksen 2003). In this paper we discuss 
possible coeval growth due to radial infall of stars and dark 
matter that both feed the black hole and establish a colli- 
sionless cusp with a self-similar distribution function. 

In this and subsequent papers (I; II; III) we will embed 
a black hole (or at least a central mass) in a distribution of 
particles that arises naturally during the formation of the 
galactic core. The 'natural evolution' is taken to be that of 
self-similarity since power-law behaviour is observed both in 
reality and in simulations. In this fashion we introduce some 
uniqueness into the form of the resulting distribution func- 
tion that describes the collisionless matter that surrounds 
the central mass. We will predict the consequent density 
cusp profile and that of the velocity dispersion variation in 
the cusp. 

In this first paper we confine ourselves to radial infall. 
The first general result on radial infall was in fact given 
in ( |Fillmore fc Goldreich, 1984[ ) and the Bertschinger (|1985p 
result follows by putting their parameter ere = 1/3. How- 
ever neither this work nor that of Bertschinger attempted to 
infer closed forms for the equilibrium distribution function. 
This was begun by Henriksen & Widrow (1995) (hereinafter 
(|HW95P 1. In (|H2006P it is pointed out that a = 9/8 (e = 3) 
yields the Bertschinger solution in its entirety, including the 



recently heralded power law of the proxy for phase space 
density. The parameter e = 3efg and it is related to the 
index called 'a' below, both by simulations and theory. 

We deduce in this paper two principal distribution func- 
tions that arise naturally. One is steady and infinite and can 
not contain a central mass exactly. However we can treat the 
central mass as a perturbation to the gravitational field and 
iterate on the Boltzmann and Poisson equations to find the 
flattest possible cusp in a steady infinite system surrounding 
a central mass. There is a special logarithmic case in this cat- 
egory that yields a Gaussian distribution function and an in- 
verse square law density cusp. Our derivation is new but the 
result agrees with the form deduced in fLyridcn-Bcll 1967| 
based on statistical mechanics. 

Our second result is most important in that it repre- 
sents a growing cusp wherein a central mass may be em- 
bedded exactly. It is the Fridmann and Polyachenko distri- 
bution function, which has not previously been identified in 
this physical context. We derive it theoretically and verify it 
with simulations. It also corresponds to a logarithmic poten- 
tial for which the self-similar index a (see below) is unity. 
The compatibility with a central mass allows us to give a 
growth rate for that mass. It is perhaps significant in light of 
the extensive orbital study of ( Van den Bosch et al. 2008) 
that we are most successful with our embedding in a time 
dependent case. For these authors suggest that the central 
mass may render the orbits chaotic and nonstationary. 

We are aware that such a radial system is unsta- 
ble to the radial orbit instability (ROI) on small scales 
l|MWH 20061 hereinafter for MacMillan et al. 2006). How- 
ever even fully cosmological simulations show an outer en- 
velope wherein the orbits are trending to be radial. More- 
over isolated halos, which also show statistical relaxation 
that is begging to be understood, show quite radial orbits in 
the envelope [MWH 2006}. In such cases the 'central mass' 
may be a central 'bulge' of dark matter that forms rapidly. 
We speculate below that the growth of this bulge in such 
an envelope may be described in terms of radial infall, and 
that this infall may continue hierarchically to smaller scales 
after relaxation by the ROI and clump-clump interactions 
(MacMillan & Henriksen 2003). Such interactions would re- 
move angular momentum from some particles in favour of 
others and so create the continuing radial infall. 

We begin the next section with the general formula- 
tion in spherical symmetry. Subsequently in section we 
discuss the various possible rigorously steady distribution 
functions (DF from now on) for radial orbits. We show 
in section 2] that the DF of Fridman and Polyachenko 
HFridmann fc Polyachenko 19"84| hereafter called FPDF) de- 
scribes a system of radial orbits that is growing self-similarly. 
This is contrasted in the same section with an infinite steady 
system for which the DF is Gaussian. After some discussion, 
we give our conclusions. 



2 DYNAMICAL EQUATIONS IN INFALL 
VARIABLES 

Following the formulation of (|H20 06) in this section we 
transform the collisionless Boltzmann and Poisson equa- 
tions to 'infall variables'. We treat a spherically symmetric 
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anisotropic system in the 'Fujiwara' form (Fujiwara 1983, 
e.g.) namely 



dr 
2 <9£ 
dr 



d®\ df 



d_ 

dr 



3 dr J dv r 

J f(r, v r ,j 2 )dv r dj 2 , 



(1) 

(2) 



where / is the phase-space mass density, $ is the 'mean' 
field gravitational potential, j 2 is the square of the specific 
angular momentum and other notation is more or less stan- 
dard. 

The 'infall variables' are a system of variables and coor- 
dinates that allows us either to readily take the self-similar 
limit or to retain a memory of previous self-similar dynam- 
ical relaxation into a true steady state. In this way we can 
remain 'close' to self-similarity just as the simulations ap- 
pear to do. These coordinates (Henriksen 2006a ; H2006) al- 
low the general expression of the Vlasov-Poisson set, but 
they also contain a parameter (a) that reflects underlying 
self-similarity. The self-similar limit is taken by assuming 
what we term 'self-similar virialisation', wherein the system 
is steady in these coordinates, although it is not absolutely 
steady since mass is accumulating in this mode. 

The transformation to infall variables has the form 
HH2006I e.g.) 



R = r e 



iT/o 



ry -2 -<A/a-2)aT 

Z = J e 



Y — v r e 

aT 



-(l/a-l)nT 



at, 



P (R, Y, Z; T) = e<V-»* T *f (r, v r , f; t) 
V(R;T) = e- 2(1/a - 1)aT -I»(r), 
Q(R;T) = p(r,t)e 2aT . 



(3) 



This transformation is inspired by the nature of self- 
similarity, which can be understood as a scaling group 
wherein each quantity scales according to its dimensions 
( |Carter fc Henriksen, 1991) 1 . The group parameter is the log- 
arithmic time T. The combinations of scaling constants 
(note that a = a/S see below) multiplying aT in the ex- 
ponential factor of each physical quantity reflect the dimen- 
sions of that quantity. When the dependence on the param- 
eter T is retained in the new variables, there is clearly no 
invariance along the scaling group motion and so no self- 
similarity. This means that the passage to the self-similar 
limit requires taking dr = when acting on the trans- 
formed variables. Thus the self-similar limit is a station- 
ary system in these variables, which is a state that we have 
described elsewhere as 'self-similar virialisation' (H W 19991 
ILe Delliou 20011 i.e. Henriksen & Widrow 1999). The virial 
ratio 2K/\W\ is a constant in this state (although greater 
than one; K is kinetic energy and W is potential), but the 
system is not steady in physical, i.e. untransformed, vari- 
ables, since infall continues. 

The single constant quantity a is the constant that 
determines the dynamical similarity, called the self-similar 
index. It is composed of two separate scalings, a in time 
and 8 in space, in the form a = a/ 5. The dimensions of 
any mechanical quantity can be expressed in the scaling 
space a = (a, 5, p), where jj, is the mass scaling. The ex- 
ponential factor of any physical quantity Q (which includes 



physical constants) is calculated as a ■ cLq where oJq de- 
scribes the quantity Q in scaling space. Thus for the velocity 
d v = (—1, 1, 0), for the distribution function df = (3, —6, 1) 
and for Newton's constant da = (—2,3,-1). In a gravita- 
tion problem the scalings a and S can express the scalings of 
all physical quantities having mass, length and time dimen- 
sions. This is because we require G to be constant under the 
scaling motion so that a- da = and hence fx = 35— 2a. This 
changes df to (1, —3) as used in equation (J3)). This procedure 
yields the constants in the exponential factors transforming 
all physical quantities in the equations (|3]). 

We assume that time, radius, velocity and density are 
measured in fiducial units r /v ,r , v and p respectively. 
The unit of the distribution function is f a and that of the 
potential is v^. We remove constants from the transformed 
equations by taking 



fo 



I 3 
; Po/V , 



4irGp r . 



(4) 



These transformations convert equations (|T|} , (J2J) to the 
respective forms 



-d T P- 

a 



--l")P+(---)^P 
a I a a 



and 



i \ i (m_ _ Z_ 

a ) + a\dR i? 3 



±JL (n 2 — 

R 2 dR V dR 



dyP-l --2] Zd z P = 



(5) 



e. 



This integro-differential system is closed by 



e = — / r,n ,iz. 



(6) 



(7) 



Until we enforce the self-similar limit (d/dT — 0) these 
equations remain completely general, because we have made 
a continuous and invertible change of variables in equation 
()3J). The merit of the transformation at this stage is only that 
it puts the expected asymptotic self-similar behaviour in the 
explicit exponential factors, while relegating the declining 
time dependence leading to this state to the transformed 
variables. These variables are strictly independent of T in 
the self-similar state. 

We will in this paper restrict ourselves to the filled loss 
cone limit of radial infall (HLcD 2002), although this is not 
the case in subsequent papers of this series. This special case 
is certainly not realistic where angular momentum becomes 
important, but it may have application on large scales and 
in the subsequent evolution in regions where angular mo- 
mentum is transformed away either by bars or other asym- 
metries. It may in any case be regarded as an introduction 
to our methods. 

To proceed we set 

P = F(R, Y; T)S(Z) = F(R, Y; T)8(j 2 ){e < - 4/a - 2)aT ) 

(SQ is the Dirac delta, not the scaling delta) which 
changes the scaling for the DF in equation ([3} to 



irf = F(R, Y, ■ T)e (1/K - 1)aT 6(f) 
while other scalings remain unchanged. 



(8) 
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The governing equations now become equation ([6| plus 
the Boltzmann equation for F(R, Y; T) in the form 



-d T F + ( - - 1 ) I- - 



Y R 



-l y 



d R F 

im_ 

adR 



Finally equation reduces to 



Q =R 2 



F dY. 



d Y F = 0. (9) 



(10) 



This completes the formalism that we will use to obtain 
the results below. 



3 STEADY CUSPS AND BULGES WITH 
RADIAL ORBITS 

In this section, we find DFs both self-similar and steady 
which are comprised of collisionless particles in radial orbits. 
We expect one mode of relaxation in collisionless cusps to be 
of the 'moderately violent' type satisfying, in terms of the 
particle energy E and mean field potential $, the relation 

dE <9$ 

This includes phase-mixing. Another mode 
QDiemand et al., 2006| ) is furnished by the presence of 
hierarchical sub-structure . The sub-structure can interact 
in clump-clump interactions that can induce relaxation on 
a coarse-grained scale ([H2 009, i.e. Henriksen 2009). 

However the temporal evolution of the system is dif- 
ficult to follow analytically even in the self-similar limit, 
so we normally look for equilibria established by the evo- 
lution. This may be either a strictly steady state in some 
appropriate coarse-grained description, or it may be a self- 
similar virialised state. 

In general one can not find a unique solution of the 
governing equations that consistently reproduce infall onto 
a central mass. We find that this is possible in one interest- 
ing case (Fridman and Polyachenko DF) of accretion onto a 
point mass that arises naturally, but not for a truly steady 
distribution around a point mass. One can allow for the 
presence of a point mass in a rigorously steady distribution 
by iterating about an equilibrium state that is determined 
initially by the central mass. This allows the central mass 
and the environment to be evolved together towards a new 
equilibrium, although normally only a single loop is feasible 
(we give two loops below as a confirmation of the continu- 
ing logarithmic behaviour). We proceed to derive these two 
results in this section. 

Using the characteristics of equation ((9]) plus the total 
derivative 



d* <9* dR^ T 

as os as 



(12) 



where ds = adT, one finds by a simple manipulation that 



(* + ») 

ds 



-2 [±-l) <"> 
a 2 a as 



isolating integral (i.e. characteristic constant), the sum of 
the last two terms must give — 2(1 /a — 1)^. This is most 
simply effected by setting d^/ds — and R8r^ = p^, 
which turns out to be a condition for both self-similarity 
and a true steady state. Here 



p = 2(1 - a). 



(14) 



so that ^ = ^oR p for some constant $ D . 

Hence, on setting £ = Y 2 /2 + ^5f, we have from equation 



ds \a 



(15) 



This variation does render E constant on characteris- 
tics (and therefore in time) as one sees by integrating to the 
form E exp (— 2(1 /a — l)s), and then by using the trans- 
formations © to find E = £ exp (2(1/ a - l)aT) = E a . An 
example of such a state is a system of massless particles 
dominated by the potential of a central point mass M», for 
which a = 3/2 and p = — 1. The similarity index a = 3/2 
reflects the presence of the Keplerian constant GM, whose 
vector (Lk = (—2, 3) and for which a ■ cLk = 0. We refer to 
this example again below. 

Equation ((9J also yields along the characteristic 



dF 
ds 



1 



-I F. 



so that with equation (|15[) 

F = F(k)\£\ 1/2 . 



(16) 



(17) 



The steady unsealed DF follows from this last equation and 
the transformations ([3]) as (a 7^ 1) 



nf = F(K)\E\ 1/2 S(f) 



(18) 



In order for this last equation to yield the energy as an 



The quantity k in equation (|18[) labels any other (be- 
sides E) possible characteristic constant, but in general 
nothing other than E is readily available. For this reason 
we discuss the DF (|18[) with F strictly constant. 

In (|HW95|) this DF was first given by assuming always 
a rigorous steady state and was shown to yield the asymp- 
totic particle distribution found by Fillmore and Goldre- 
ich ( |Fillmore~&: Goldrcic h, 1984[ | .Their result followed from 
a direct integration of particle orbits in radial infall. The 
example given in l|HW95[l for comparison purposes corre- 
sponds to the choice of our current index a = 18/17. In 
general, a — 3e/2(e + 1), where the initial density profile 
is oc r _E . The initial system is infinite for e < 3, which is 
a < 9/8. In l|HW 1999[) this DF was argued to be the natu- 
ral state for steady self-similarity. Thus the DF (|18[) with F 
constant probably describes a steady system of self-similar 
radial orbits characterized by the index a. For this reason 
together with the numerical evidence discussed presently, we 
discuss the implications of this DF here in more detail, bear- 
ing in mind the possibility of embedding a central black hole 
or other spherical mass. In (|HW 1999|) weak evidence was 
presented to show that the IF] 1 / 2 law did appear near the 
end of the infall for the most tightly bound particles with 
negative energy. These would be closest to being described 
as occupying a rigorous steady state. A similar result was 
found by MacMillan (jMacMillan 2006j) for the most tightly 
bound particles even while infall continued. The case was 
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Figure 1. A shell code evaluation of the DF (Le Delliou 2001), evolved from a system with initial density p oc r — 15 . The first fit is with 
a cut off power law (F oc E~ v e~ E / Ea with p ~ 1/2, E c ~ — 10~~ 2 ; upper left panel), while the second fit is just a power law (lower left), 
confirming (HW 1999). The cut off is confirmed by higher resolution in the DF (upper right) and in density of states g(E) (lower right). 



reinforced by additional calculation of initially infinite sys- 
tems by Le Delliou (ILe Delliou 2001)1 (see Fig.Q]). The figure 
shows another example of the Fillmore & Goldreich (1984) 
problem, wherein the DF (JTHJ) is a reasonable fit over most 
of the energy range. The cut-off is probably numerical in 
this case since numerically the system is ultimately finite. 

However the DF (|18p does not fit the complete energy 
distribution found in high resolution simulations of radial or- 
bit growing isolated halos (MacMillan 2006, e.g.) in a state 
of self-similar virialisation (|HW 199 9). In such a state, in- 
fall continues. We reserve the explanation of this behaviour 
to the next section where it involves the exceptional case 
wherein a = 1. 

Continuing for the moment with the discussion of equa- 
tion (|18[) we observe that an upper energy cut-off is required 
for finiteness in the positive energy case (a < 1 when the po- 
tential increases outwards), while the cut-off is zero in the 
negative energy case (a > 1 when the potential decreases 
outwards). Moreover we note that one can add an arbitrary 
constant E to E in this DF, which reflects the arbitrary 
constant in the potential. For negative energy the DF would 
appear as F(E - E) 1/2 for E < E a < 0. This was the 
type of fit used by MacMillan (2006) to fit his simulations. 
The DF decreases to zero at E and increases to negative 



energy. A large positive constant E and E < E so that 
/ oc (E — E) 1 ^ 2 would express a positive energy cut-off at 
the value E and the DF would increase towards zero energy. 

The potential and density pair for these rigorously 
steady models take the forrrQ (a 7^ 1) 

* = *o-R (2_2a) , e = 2(3-2a)(l-a)* ir 2a . (19) 

In a very recent submission Amorisco and Evans 
l|Amor iscofeEvans 2010) prove that the power-law relation 
between the galactic half-light radius and the central ve- 
locity dispersion in dwarf ellipticals requires a power law 
potential to be valid. This arises naturally here. 

In the example of a dominant central point mass, insert- 
ing the index a = 3/2 in the above pair yields a point mass 
surrounded by massless particles. The massless particles 
may be distributed in any manner, but since djv = (0, —3) 
dimensional analysis requires that the number density N 



1 We have used our current notation when using the results of 
previous papers, in which 8,X,S in (HW 1999) are respectively 
1 /a, R and SR 2 . In AHW95I I 8 there becomes 1/(1 — a) in current 
notation. 
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should vary as N oc N (R)e~ 3ST = N a (R)R 3 /r 3 . In a rigor- 
ously steady state this should be time independent so that 
N (R) oc R~ 3 and hence N oc r~ 3 . Such a halo could exist 
outside any dominant mass as was discussed in (HW 1999). 
Thus it might surround a central point mass or indeed be 
the diffuse halo around a bulge containing most of the sys- 
tem mass. However this is only a limiting behaviour and 
does not include the transition region. This region interests 
us particularly in the context of central black holes. 

The direct density integral over the DF (|18[) for negative 
energies yields for p 



(20) 



_ txF |$| 

Since this is linear in the potential, one can readily include a 
central mass by iteration. We may begin with a point mass 
potential for $ in the density (|20[). and then use the Poisson 
equation to obtain a new potential in a form that is no longer 
self-similar. This yields 

M* + C 2 (l + lnr) 



$1 



(21) 



whence follows a new density by (|2U0 . The constant M* 
would be the mass inside lnr = — 1 (regarded as a point 

mass) while C 2 = (itF/^/2^ A/*. There is only a logarith- 
mic modification to the r~ 3 law at large r where the iteration 
should apply. In effect the iteration yields a singular pertur- 
bation series at r — because of the diverging potential and 
hence energy. Therefore we arbitrarily cut off the series at 
small r and 'renormalize' the central mass to the mass inside 
this cut-off radius. The next loop of the iteration gives 



<3> 2 = 



M.(l- 



2M- 



) 



C 2 (l + 



2C 2 ,(l + hir) 



hr 



2M, 



where again we have 'renormalised' at lnr = —1. Putting 
this back into equation (|20[) yields p 2 as r~ 3 flattened only 
by logarithmic terms at large r (|HW 199 9)) as expected. The 
large scale r~ 3 density profile does not fit the bulge simula- 
tions inside the NFW ( |Navarro, Frenk fc White 1 996 ) scale 
length, but it does describe the halo region outside a central 
bulge of mass M* l|HW 19990 . 

Since the density is linear in the potential we may also 
solve for a self-consistent cusp having the DF 1)180 by let- 
ting the potential be determined by the Poisson equation. 
Working in transformed variables we find 



* = -AR p - 



BR P 



where A, B are arbitrary real constants > and 



P± 



V2 



(22) 



(23) 



By letting F — > we see that p_ is the power that 
should be taken near the centre if we wish to create a strong 
central mass concentration. It tends to —1 in this limit while 
p+ tends to zero. Hence we set B = in this limiting domain. 
The potential then satisfies our basic condition ()140 with 
a new self-similar index. This is given by a = 1 — P-/2 
according to equation (|14[) . that is explicitly 



- + i. /- —F 

4 + 2\/4 ^2 



(24) 



Equation (|19[) now gives the inner cusp density law as 



e = \p. 



(-2+P-) 



(25) 



This can not be flatter than R ' , which appears only for 
the ' maximum bulge' for which F = 1/(2^^2). At large r 
the term in p + dominates, and the behaviour tends to r~ 2 
for F small. 

In the context of dark matter simulations such a steady 
halo of radial orbits could describe the region just beyond 
the NFW scale radius (which we take to form the 'bulge'), 
based on the density profile alone. It is not stable in a 
strictly steady state according to the usual Antonov criteria 
(Binney & Tremaine 1987, e.g.), unless the energy is nega- 
tive. Consequently we do not expect it in central regions 
where a < 1 and the energy is positive (with a central zero: 
the potential increases outward according to Eq. 1190 . The 
radial velocity dispersion is v 2 = |<I ) |/2. 

This concludes our study of the general, steady, spher- 
ically symmetric, DF for power law distributions of radial 
orbits. The main justification for the study is that it per- 
mits definite conclusions. The rigorously steady DF (1181) is 
not permitted to contain a black-hole or other centralised 
spherical mass without perturbation. By iteration we show 
that it can only yield logarithmic corrections to an inverse 
cube law. This suggests that a black hole will engender time 
dependence in its surroundings. 

The DF (|180 also permits a self-consistent exact solu- 
tion for the potential and density in the cusp (see e.g. equa- 
tion [25}. However it fails to produce a sufficiently flat cusp 
of dark matter or stars to explain observations of the Milky 
Way cusp of stars. At large distance it yields (at flattest) an 
inverse square law density. 

For these reasons we turn in the next section to a time 
dependent case which is much more promising. It allows us 
to calculate the growth of a central mass due to infall of 
collisionless matter on radial orbits. 



4 THE LOGARITHMIC CASE 

The case a = 1 is obviously special and, as it turns out, 
rather important. We return to the equation (|130 and ob- 
serve that we also have an integral in the steady state if 
ip — ipo In R. For in that case the equation integrates to 
£ + ip aT = k, where n is constant on the characteristic. 
Equation ()160 then requires that in general F = F(k), so F 
is also constant on the characteristic. With a — 1 we note 
from equation ([3| that = $, Y = v r and so £ = E. How- 
ever k = vl /2 + Vo In R + aT = v 2 r /2 + ip In r = E and 
moreover irf = F. Hence our conclusion is for the moment 
only that 



7T/ = F(E), 



(26) 



in this case. 

By taking the potential to be logarithmic, ^ = ^ In R, 
we require by equation ([6]) that & — ^ /R 2 and hence, by 
the appropriate members of the set ([3]), that p = ^> /r 2 . But 
for consistency we must also have 



1 



F{E) dv r 



(27) 
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We must therefore find a DF F(E) which satisfies this equa- 
tion. In effect, the integral over the particle velocities must 
be a constant independent of the logarithmic potential. 

To find such a DF we convert the integral to an integral 
over energy in the normal fashion and write our consistency 
condition as 

F(E) 



V2 



dE = * . 



(28) 



We might expect a power law form for F(E) on general 
grounds, and given this a brief experimentation shows that 
the most general form for F(E) in this case may be written 
as (we suppose negative energy to ensure convergence and 
E < E < where E is an arbitrary constant energy) 

F(E} =7m=w> {m 

This may also be inferred as the unique solution with a 
finite energy range by recognizing that equation (|28|) is 
really a simple form of Abel's integral equation (see e.g. 
(Binncy & Trcrnainc 1987) in the first or finite form, p65f) 
whose solution is equation (|29|) with 

K = 7t- m 

This may be checked by direct evaluation of the integral to 
find that p — fy /r 2 as required. 

This result holds only where E < and hence where r < 
r , where r is an arbitrary scale. However since we have not 
actually set a fixed scale in the problem (which would entail 
setting 8 = 0), we have implicitly implied that r/r = R/R . 
Thus we have assumed that r — R e aT = aR t. We may 
take R constant (dimensionless) so that there is a residual 
time dependence because the outer boundary expands ac- 
cording to r = R e aT . This implies that the mass inside R 
and indeed inside any fixed R is growing as M = 4n^ Rat. 

We have thus succeeded according to the above 
in deriving the Fridmann and Polyachenko DF 
( |Fridmann fc Po lyachenko 1984) as the unique result of 
time-dependent radial accretion of a growing inverse-square 
density 'bulge'. This is one of our major conclusions. 

A numerical measure of the DF in radial self-similar 
continuing infall was made by MacMillan (MacMillan 2006). 
Instead of the steady DF, the DF of Fridman and Poly- 
achenko (Fridmann & Polyachenko 1984) is found to predict 
accurately all of the measured quantities as in the accompa- 
nying figures. These include an inverse square density law 
and a power law pseudo-phase-space density of w —1.5, 
but there are logarithmic corrections to the power law as 
can be calculated from the FP distribution function. The 
pseudo-phase-space density power is flatter (MWH 2006] i.e. 
MacMillan, Widrow & Henriksen 2006) than is generally 
found in full cosmological simulations. These results are il- 
lustrated in figures and ([3]). 

This DF (Fridmann & Polyachenko 1984) used to make 
the fits in figures © and @ is 

f = (-E + Eo y^ ^ 

for E < E < 0, and r < r$ (where "l?(r/) = E a ) and 
zero otherwise. This is just as we inferred above. The den- 
sity profile is r~ 2 and the potential is logarithmic. The log- 
arithmic variation of the velocity dispersion together with 




Figure 2. We show the Fridman and Polyachenko fit to the 
mass distribution dM/dE = f(E).g(E), density of states g(E), 
and the phase space distribution function f(E). The figure is 
based on the radial simulations of an isolated dark matter halo 
by (MacMillan 2006). The system is maintained in self-similar 
virialisation by steady accretion. The fits use equation J 3 IP with 
K = — _E /(4v / 2~7r 3 ) and E Kb —80 in machine units. 



the inverse square density profile accounts for the pseudo 
density approximate power law found in the simulations 
(|MacMillan 2006|) . 

The persistence of this DF is undoubtedly due to the 
strict proscription of non-radial forces in the simulations. It 
is not linearly stable by the Antonov criteria for E < 0. 
When this proscription is relaxed (MacMillan 2006) shows 
that the equilibrium Fridman and Polyachenko DF is sub- 
ject to the radial orbit instability. It may require continual 
non-equilibrium excitation as provided by steady infall to 
be realized. 

The unique feature of the distribution function (|31|) is 
that the density is independent of the potential. Hence one 
can simply add a point mass potential to the logarithmic 
bulge value and the density will remain ^ r~ 2 . In the case 
of a true absorbing central black hole one should only permit 
negative radial velocities in the system. This means that K 
in equation (|30p should be multiplied by a factor 2. The 
velocity dispersion however goes as v}. = $ — E \ and we 
may take <f> = -M*/r + V2nK In r + E a . MacMillan (2006) 
finds a good fit to this radial dispersion in his simulations, 
but without a central mass. 

The actual growth of the black hole mass will be sim- 
ply that of the general self-similar mass growth as discussed 
above. That is 



M,(t) = M„(0)e° 



M t (0)at. 



(32) 



Its radius will be growing according to the same law. 

One can not however take seriously the growth of a 
black hole due to radially infalling material from cosmologi- 
cal distances, since the material there can not know the ac- 
tual location of the black hole and the radial infall is subject 
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Figure 3. We show the the Fridman and Polyachenko fit to the 
mass density, the velocity dispersion and the 'pseudo-phase-space 
density' for the same simulations by MacMillan. 



to the radial orbit instability. However, as we will speculate 
in the discussion section, such growth may apply to a hi- 
erarchy of 'central' masses extending to ever smaller scales. 
That is M* might be successively the bulge of a galaxy, the 
core, the nucleus, and so on to the black hole. In each case 
there must be a way of scattering orbits into the essentially 
radial loss cone. Such scattering may be due in part to the 
formation of a bar by the radial orbit instability itself, or to 
clump-clump interaction. 

It is worth contrasting the above infall time-dependent 
DF with the rigorously steady system of radial orbits in 
spherical symmetry. This case was inadvertently omitted in 
the discussion by (HW95) and so we pause to present it in 
our general style as another example of the method, similar 
to that above but in a strictly steady state. 

The main difference with the time dependent case is 
that that the scaling motion must be taken in space rather 
than in time. This means that we use the variable R, where 

a SR 



Sr so that dR/dr 
f = 



e SR . In addition we write 
F8(ve)8(v^, 



(33) 



where F satisfies [integrate the general steady CBE in spher- 
ical symmetry (Binney & Tremaine 1987) over vg and v$\ 



OF 
dr 



9$ dF_ + 2Fv r 
dr dv r r 



= 0. 



In addition we have the Poisson equation as 



_d_ 

dr 



,d$ 



= p= I F dv 



(34) 



(35) 



The dimension vectors in our usual scaling space are 
d F = (1,-4,1), d v = (-1,1,0), dp = (0,-3,1) and d* = 
(—2,2,0). Recalling that /i = 35 — 2a the vector for F be- 
comes (—1,-1), the vector for p becomes (—2,0) while the 
others remain unchanged in the reduced (a, 5) space. How- 
ever we will only consider the case a = 1 or a — 8 since the 



other cases were discussed in \ 1 1\\ and reduce to the DF 
(|18[) . When a — 8 the dimension vectors reduce in delta di- 
mension space to (—2), (—2), (0) and (0) respectively. Hence 
the equivalent of equation ((3]) is 

SR 



dr : 



e~ 2SR P (R, Y) =F{r,v r 



Y 



(36) 



*{R) = $(r,t), 



e~ 2SR Q {R) = p{r). 

If self-similarity is enforced normally, then in this case 
8r — when it acts on the scaled variables. However this 
can not apply to '3/ (7?) = <& since then equation Q34p becomes 
trivial. The answer lies in the realization that in this case 
the potential is logarithmic in r, rather than being a power 
law. The condition a — S requires that there be a constant 
*l/ with dimensions of the potential per unit mass, just as 
in the time dependent case. These conditions are satisfied 
here by setting 

V(R) = ^o(SR) =e * hi(5r = $( r ). 

Since only 8r^ appears in the problem, the self-similar re- 
quirement of independence of R is maintained for F. 

A direct substitution of the transformation (|36p plus 
the form of the potential into equation (|34p reduces it to 

d$ 

Yd R P ~ -TB d Y P = 
dR 

and hence after solving by characteristics 
P = P(E) 

Consequently we find finally from F = Pe~ 2SR = P/(Sr) 2 



the steady density in the form 



where 



1 

8 2 r 2 



E 



P(E) dv r 



(37) 



From the Poisson equation l|35[) we obtain 



and this must agree with the integral over F. Letting both 
inward and outward going particles be present one finds that 

V2 



P(E) 



dE. 



vW-S) 

The problem is now reduced to the same situation that we 
had in the time dependent case found earlier (|27l and follow- 
ing). In order to have a consistent (with the Poisson equa- 
tion) inverse square density law, only the steady Fridmann 
and Polyachenko DF may be invoked, and that only for nega- 
tive energies wherein $ < E < E . Thus we derive the finite 
and steady Fridman and Polyachenko result as quoted for 
example in ( |Binney fc Tremaine 1987). 



5 DISCUSSION AND CONCLUSIONS 

In this paper we have presented unique distribution func- 
tions based on the self-similar evolution of a gravitating 
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system of radial orbits. The rigorously steady state corre- 
sponds to an infinite system while the time dependent sys- 
tem is finite and growing. Both of these systems have been 
confirmed by numerical simulations, with the confirmation 
of the growing mode being most dramatic. 

Both systems have been derived coherently from the 
basic equations by using a convenient formulation of self- 
similarity. We have considered how these distributions are 
affected by the presence of a central spherical mass, which 
may not necessarily be point-like. The steady system is per- 
turbed outside the central mass to an inverse cube law with 
logarithmic corrections as we discuss further below. The 
growing time-dependent case forms a perfect inverse-square 
density cusp that can contain a central mass and is grow- 
ing proportionally to the time since formation. The material 
with which we are concerned is collisionless and so may refer 
to dark matter in an early stage of formation or to stars at 
a later stage. The restriction to radial orbits in this paper is 
removed in subsequent papers, although both in theory and 
in the numerical simulations a bias to radial orbits remains 
on the larger scales. 

For this latter reason we think that the results of this 
paper may be relevant in reality for two reasons. In the first 
instance the distribution functions may develop outside a 
rapidly formed large central bulge wherein angular momen- 
tum has been important. This may correspond to the NFW 
core radius, and our detailed results are compatible with the 
density profile of the simulations outside this radius. The 
second plausible scenario is that, on much smaller scales 
surrounding an actual black hole, various instabilities may 
allow radial infall of stars and dark matter onto the cen- 
tre. This fits the time-dependent infall model and allows the 
black hole to grow proportionally to the time from the onset 
of the radial infall. We detail these results in what follows. 

In section Q on the steady state we deduced the steady, 
spherically symmetric DF with radial orbits that yields infi- 
nite systems with power law profiles. This was found pre- 
viously but we rederive it here as a limit of the time- 
dependent equations. We have perturbed this DF by embed- 
ding a central mass. The DF is universal for these systems 
(/ = K\E — E\ x ^ 2 ) but the potential-density pair depends 
on the self-similar index a which in turn depends on domi- 
nant constants or boundary conditions. If it is a memory of 
a cosmological fluctuation profile, then a — 3e/(2(e + 1)). 

This steady DF can not contain a central mass with- 
out being perturbed, except for the Keplerian limit wherein 
a — 3/2 and there is a mass- less cusp with an inverse cube 
density profile. Iterating the equations to find the pertur- 
bation produced by an embedded central mass predicts a 
transition region in the halo of the mass that is an r -3 pro- 
file modified only by logarithmic corrections. It does not 
correspond to the power law density of stars found close to 
the black hole in the galactic centre ( [Gillessen et al., 2 009), 
but may apply outside a bulge mass confined to the NFW 
scale radius. 

Since the density is linear in the potential we were 
also able to find a solution that broke the self-similarity by 
seeking non-power-law solutions of Poisson's equation (pure 
power laws only exist in the limits) with the density (120[) . 
The density profile of the cusp can not be flatter than r -2 ' 5 
near a central point mass, but it may be as flat as r~ 2 at 
large scales in the low density limit. Such behaviour is too 



steep to explain the cusp observed around the Sagittarius A* 
black hole. All possible descriptions of a steady system of ra- 
dial orbits are clearly excluded by this observation. These 
descriptions are not excluded outside a central bulge how- 
ever. 

In the section concerning the logarithmic potential ex- 
ception we showed that it corresponded to continuing time- 
dependent infall. We found both theoretically and numeri- 
cally that the Fridmann and Polyachenko DF is the distri- 
bution established by continuing radial infall. Remarkably, 
it allows a growing central point mass (or indeed a bulge 
mass on a larger scale) to be embedded in the infall self- 
consistently. In both the steady and the time-dependent 
cases the density profile is r~ 2 . The time-dependent log- 
arithmic case allows a consistent calculation of the cen- 
tral mass growth rate according to M» oc t. Such a 
growth rate from a surrounding envelope was found in 
(MacMillan & Henriksen 2003). This simulation used a true 
N-body code that treated particle-particle scattering, which 
scattering led eventually to radial infall of some of the par- 
ticles. 

The growth rate of a central mass (i.e. a collisionless 
concentration, not a true black hole) from a reservoir of ra- 
dial orbits is zero if there is a rigorous steady state. The 
growth rate is not zero if the central mass is a black hole, 
since then the outward bound radial orbits are suppressed. 
In that case however the steady state is only an approxi- 
mation except in an infinite system. For a finite system the 
true timescale would be the free-fall time of the bulk of the 
accreting mass. 

One is inclined not to take these growth estimates seri- 
ously, since they simply assume an endless supply of radial 
particles, and hence an arbitrarily filled loss cone. In fact 
the radial alignment required to hit a growing black hole 
from a few hundred parsecs is at least one part in 10 8 to 
10 10 depending on the mass of the black hole! This sug- 
gests that instead (MacMillan & Henriksen 2003, see e.g.) 
the actual growth involving radial orbits may be by way 
of a multi-stage process. In the first stage, radial orbits 
accrete from the galactic halo to form a bound spherical 
bulge of intermediate size, due to finite angular momentum 
about the centre (MacMillan & Henriksen 2003). They are 
trapped there either by the usual mechanism of self-similar 
infall as the potential increases in time with increasing inter- 
nal mass, or by dissipative interactions. If there is substruc- 
ture in the collisionless matter (e.g. stars and dark mat- 
ter clumps), then these are able to produce dissipational 
collisions. Ultimately these collisions and tidal interactions 
can lead to a more gradual growth of a more central mass 
(MacMillan & Henriksen 2003, e.g., in the Carnegie meet- 
ing)- 

Recently a high resolution study (|Stadel et al 2009)) of 
sub-structure in an isolated halo revealed that this sub- 
structure disappears in the inner few parsecs. this coincides 
with the region where the halo is becoming spherical and 
where the density power-law is flattening to less that 1. The 
interactions leading to the sub-structure disappearance may 
well lead to relaxation. 

In addition the radial orbit instability can lead to the 
development of a bar l|MWH 2006p . This bar can then trans- 
port angular momentum away from the bulge by the ejection 
of particles. Such 'interrupted accretion' may repeat several 
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times on the way to the actual central object. The rapid 
accretion of a bulge is in fact the way in which dark matter 
halos are thought to grow UZhao et al., 2003|lLu~et~aI~2 006) 
initially. This is then followed by a slower growth phase. The 
DF (|31[) can be used to describe the environment of the cen- 
tral mass on each scale of the interrupted cascade. 

In this connection we refer to the work of Mutka 
(Mutka 2009) on gravitationally lensed galaxies with double 
images. He concludes that there are two classes of density 
cusps with the larger sample (about 80%) showing a loga- 
rithmic density slope of w —1.95 well inside the NFW scale 
radius. The other 20% show this slope as « —1.45. These 
may be unresolved triple image lens and ,if so, the measured 
value should be rejected. 

Mutka's result is a measure of the total mass distribu- 
tion rather than just the dark matter. Perhaps we are seeing 
enhanced relaxation in the mixture of stars and dark mat- 
ter, that leads towards an isothermal cusp, rather than the 
shallower cusps of the dark matter simulations. It is signifi- 
cant that this inverse square slope is also frequently found by 
direct dynamical modeling of galaxies (|van der Marel 2009 ). 

However an inverse square slope is not restricted to a 
system of purely radial orbits as the isotropic isothermal dis- 
tribution shows. In a subsequent paper we survey anisotropic 
distribution functions in spherical spatial symmetry that 
also have a self-similar memory. Some of these also provide 
an inverse square density profile. 

The significance of the hierarchy of co-evolving struc- 
tures is that there will always be a mass correlation between 
them. Thus if the mass concentration derives its ultimate 
mass Mm from a halo of radius rn, while r s encloses the 
mass that forms the ultimate bulge mass M s then 



This assumes the pure inverse square density law, which 
might in fact have a logarithmic correction. In the subse- 
quent paper, we shall find a slightly more general correla- 
tion that involves the self-similar memory. Taken at face 
value this simple relation gives rh/r s ~ 100. 

This paper comprises a systematic derivation of the dis- 
tribution functions that arise through self-similar evolution 
of gravitating systems of radial orbits. In addition we have 
studied the effects of a central mass on these distributions. 
We do not find density profile as flat as those proposed 
in (|Merritt fc Szell 20 06 ) and (|Nakano fc Makino 1999]) as 
a result of density scouring by merging black holes. Thus 
our calculations probably do not describe the cusp around 
Sagittarius A*, although this does not exclude other applica- 
tions. Subsequent papers in this series address this problem 
by allowing anisotropy in velocity space. 
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